This notebook will run and then plot using Delay and Sum photoacoustic image reconstruction method.
The Delay and Sum algorithm works by synchronizing the arrival times for an array of transducer elements, summing their values, and averaging. The result represents the initial pressure value for a discrete point in space. Repeated over a region in space and one can create an initial pressure image from the photoacoustic energy transfer function.
Import Packages
from scipy.io import loadmat
import numpy as np
import plotly.express as plt
import time as time
import math as mt
from numba import cuda
Kernel Function: InitIndexMat
Inputs:
xImageLocs: all x-values for each reconstructed pixel
yImageLocs: all y-values for each reconstructed pixel
sensorLocs: [x, y] position for each element
indexMat: empty matrix for each the corresponding index of each element's pixel
signalLength: length of rfdata's first dimension
sos: medium speed of sound
dt: inverse of the sampling frequency for rfdata
Output:
N/A
Purpose: To calculate the index matrix for each pixel and element. It assumes that indexMat will be used to index a flattened signal matrix
@cuda.jit
def InitIndexMat(xImageLocs, yImageLocs, sensXloc, sensYloc, indexMat, signalLength, conversion):
layer, locx, locy = cuda.grid(3)
if layer < indexMat.shape[0] and locx < indexMat.shape[1] and locy < indexMat.shape[2]:
indexMat[layer, locx, locy] = (mt.sqrt((xImageLocs[locx, locy]-sensXloc[0,layer])**2+
(yImageLocs[locx, locy]-sensYloc[0,layer])**2)/conversion
+ (signalLength*layer))
cuda.syncthreads()
Class Name: VariableSetup
Inputs:
filepath: exact file location for the .mat datafile
xlims: min and max values for the x-axis in the reconstruction region
ylims: min and max values for the y-axis in the reconstruction region
res: pixel resolution in the reconstruction region
Attributes(memory host):
.sens_count(CPU): number of transducer elements
.rfdataALL(CPU): If more than one frame is present this holds all frames
.rfdata(GPU): time vector with pressure values for the first frame recorded by each element
.sensLocs(CPU): [x, y] position for each element
.c0(CPU): medium speed of sound
.dt(CPU): inverse of the sampling frequency
.xMat(CPU): all x-values for each reconstructed pixel
.yMat(CPU): all y-values for each reconstructed pixel
.allIndex(GPU): indeces corresponding to each pixel in the layered delay matrix for all elements
l x m x n, l - corresponding to each element layer, m x n - corresponding to the reconstruction region
Methods:
__init__(self, filepath, xlims, ylims, res):
initialize a VariableSetup object with attributes above and one conditional;
For multiple frame rfdata, only the first frame is sent to the GPU.
SetupGPU(self, tb_shape):
create the index matrix for all image reconstruction frames, calls above "InitIndexMat" kernel function on cuda.
Inputs:
self - VariableSetup object
tb_shape - threads per block shape for cuda kernels
Return:
blocks_grid - shape of the blocks per grid for future image reconstruction
Purpose: To hold the necessary information for calculating DAS images given a constant system setup and speed of sound.
class VariableSetup:
def __init__(self, filepath, xlims, ylims, pixLength):
dat = loadmat(filepath)
self.sens_count = np.squeeze(dat['ele'].astype(int))
if np.array(dat['rfdata']).ndim > 2:
self.rfdataALL = np.array(dat['rfdata'])
self.rfdata = cuda.to_device(self.rfdataALL[:,:,0])
else:
self.rfdata = cuda.to_device(np.array(dat['rfdata']))
self.sensXlocs = np.array(dat['x0']) * 1e-3
self.sensYlocs = np.array(dat['z0']) * 1e-3
self.c0 = dat['sos'] * 1e3
self.dt = 1 / (dat['fs'] * 1e6)
self.xMat = np.transpose(np.tile(np.linspace(xlims[0], xlims[1], pixLength), (pixLength, 1)))
self.yMat = np.tile(np.linspace(ylims[0], ylims[1], pixLength), (pixLength, 1))
self.allIndex = cuda.device_array((int(self.sens_count),pixLength,pixLength),dtype = 'uint64')
def SetupGPU(self, tb_shape):
threads_block = tb_shape
blocks_gridx = mt.ceil(self.allIndex.shape[0]/threads_block[0])
blocks_gridy = mt.ceil(self.allIndex.shape[1]/threads_block[1])
blocks_gridz = mt.ceil(self.allIndex.shape[2]/threads_block[2])
blocks_grid = ( blocks_gridx, blocks_gridy,blocks_gridz)
InitIndexMat[blocks_grid, threads_block](self.xMat,self.yMat,self.sensXlocs,self.sensYlocs,
self.allIndex,int(self.rfdata.shape[0]),float(self.c0*self.dt))
return blocks_grid
Kernel Function: FrameImaging
Inputs:
indexArr: Array from VariableSetup class .allIndex
valueArr: Flattened array from VariableSetup class .rfdata, only one "frame" at a time.
outputArr: Array with the same shape and size as the indexArr, used to store initial pressure values for each pixel
Output:
N/A
Purpose: To calculate the initial pressure values for each pixel averaged by each transducer element.
@cuda.jit
def FrameImaging(indexArr, valueArr, ImgArr):
layer, locx, locy = cuda.grid(3)
if layer < indexArr.shape[0] and locx < indexArr.shape[1] and locy < indexArr.shape[2]:
ImgArr[layer, locx, locy] = valueArr[indexArr[layer, locx, locy]]/indexArr.shape[0]
cuda.syncthreads()
xlimits = [-0.015, 0.015]
ylimits = [-0.015, 0.015]
pixNum = 960
file = "D:\Documents\PSU\Research\DelayAndSum\Circ_Array\Leaf_exp\leaf.mat"
DASvars = VariableSetup(file,xlimits,ylimits,pixNum)
tpb = (4,16,16) #Unique to GPU, mine fits this shape well
bpg = DASvars.SetupGPU(tpb)
imageGPU = cuda.device_array((int(DASvars.sens_count),pixNum,pixNum), dtype = float)
C:\Users\paul\anaconda3\lib\site-packages\numba\cuda\cudadrv\devicearray.py:790: NumbaPerformanceWarning: Host array used in CUDA kernel will incur copy overhead to/from device.
warn(NumbaPerformanceWarning(msg))
Create the image- This can be run in a loop with multiple frame rfdata, re-assigning rfdata for each frame
GPUstart = time.time()
FrameImaging[bpg, tpb](DASvars.allIndex, DASvars.rfdata.ravel(order='F'), imageGPU)
onGPUend = time.time()
image = imageGPU.copy_to_host()
offGPUend = time.time()
print('Time taken on GPU: ', onGPUend-GPUstart, ' [s]')
print('Time taken total: ', offGPUend-GPUstart, ' [s]')
Time taken on GPU: 0.20968008041381836 [s] Time taken total: 1.0188746452331543 [s]
Show the Reconstruction- optional to view the reconstructed image
fig = plt.imshow(np.mean(image,axis=0),
labels=dict(x='mm', y='mm'), x=(DASvars.xMat[:, 0]).astype(str), y=(DASvars.yMat[0, :]).astype(str))
fig.show()